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The activity of cortical neurons is correlated, and this correlation af- 
fects how the brain processes information. What are the neural circuit 
mechanisms underlying correlations? Shared synaptic input is necessary 
O^i to synchronize neurons, but it may not be sufficient if neural circuits are 

heterogeneous. We study the neural circuit mechanisms of correlations 
cn ■ by analyzing a network model characterized by strong and heterogeneous 

interactions. We consider strong excitatory and inhibitory connections: 
the excitatory input drives the fluctuations of neural activity, which are 
counterbalanced by the inhibitory feedback. While the strong excitatory 
input tends to correlate neurons, the inhibitory feedback strongly reduces 
the correlations. Crucially, the heterogeneity of synaptic connections is 
necessary for the inhibition of correlations. Under these conditions, we 
found that correlations are positive and of magnitude K~2^ where K is 
the number of connections to a neuron. The result agrees with anatomical 
and physiological measurements in the cortex: a cortical neuron receives 
about K ~ 10^ connections, and the measured correlations are about 10~^. 



Introduction 

Simultaneous measurements of the activity of multiple neurons have shown sig- 
nificant correlations, and this observation has stimulated the debate on whether 
and how correlations contribute to neural computation. In principle, correlations 
allow robust signal processing, because redundancies across neurons can be ex- 
ploited to separate the signal from the noise [H |3T] . Experimental studies of the 
cerebral cortex suggest that correlations improve decoding of stimuli [18] , but it 
remains unclear whether a parsimonious decoder should rely on correlations [2]. A 
challenge to this hypothesis is the observation that correlations are reduced when 
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animal subjects are actively engaged in discrimination [9l[T0], and even when they 
simply start moving [32]. In addition, neurons with similar responses to stimuli 
show higher correlations [44l|26l[27l|4l[12l|2l|36l[23l[39ll22l[l5l[S], implying that 
discrimination of stimuli is worsened by correlations [H [311 SOI SSI El [IE] • Another 
caveat is that the neural code is largely unknown, and if the "noise" measured 
in physiological studies encodes some signal then any correlation would decrease 
the available information [29] . 

Besides the possible function of correlations in signal and information process- 
ing, their physiological causes remain unclear. It has been shown that the corre- 
lation between nearby neurons is driven by their shared synaptic input [25| 132] . 
However, a quantitative understanding of the circuit mechanisms regulating cor- 
relations between cortical cells is still missing, and the goal of this study is to 
determine the dependence of correlations on different properties of the neural 
circuitry. The measured correlation between neurons depends on different fac- 
tors and varies across studies [TT]: it increases with the proximity of neurons 
[27] [T2| [39] [T5| |23], their activity [13] and the timescale on which the activity 
is integrated [ll|3l[T2l|2l|23l[3Sl|22l[28]. Here we do not consider the effects 
of distance and firing rate, and we study the correlation of the instantaneous 
activity, calculated on the timescale of a single action potential (about 1ms). 
Fig.l shows the correlation measured in eight different studies as a function of 
timescale; When measured at short timescales, the correlation is positive and of 
the order of 10~^. 

Previous modeling studies of neural circuits have found that the mean corre- 
lation between neurons is small, of the order of A^~^, where A^ is the number of 
neurons in the network. Small correlations have been observed, not surprisingly, 
in networks characterized by weak connection strengths [T71 [S]. More surpris- 
ingly, the same result has been obtained in the case of strong connections, the 
high- conductance state [HIST], provided that the network includes a strong in- 
hibitory feedback [351 [21] • Here we provide an analytical study of correlations, 
and we confirm both the observed small correlation, and the crucial effect of the 
inhibitory feedback in reducing it; However, we find that correlations are of mag- 
nitude K~^, where K is number of connections received by a neuron, instead of 
the A^~^ observed in previous studies. We argue that our result is consistent with 
empirical observations, since a cortical neuron receives about K ^ 10^ connec- 
tions and the correlation is ~ 10~^. In addition, we show that this result crucially 
depends on the heterogeneity of connection strengths. 

Methods 

The model consists of a neural circuit of A^ neurons, receiving input from Next 
external neurons. The dynamics of neuronal activity is described by a linear 
firing rate model: each neuron integrates the signal from other neurons, weighted 
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Figure 1: Mean correlation across neuron pairs plotted vs the timescale on 
which correlations are measured. Re-plotted from eight experimental studies of 
the cortex (legend). Two studies provided not only the mean correlation but also 
the Standard Deviation (SD, inset). 



by the synaptic connection strength. The dynamics of the model is described by 
the equation 

r^ = -x.{t) + f: G.,x,{t) + Y: Gfxfit) (1) 

i=i i=i 

where Xi is the activity of neuron i in the local circuit, and Gij is the strength of 
the synaptic connection from neuron j to neuron i. The external (feed- forward) 
input to the circuit is provided by the activities x^^*, and the synaptic connection 
from the j-th external neuron to the i-th local neuron is given by the strength 
GfJ*. All neuronal activities evolve in time, while the connectivity matrices G 
and Gext are fixed. 

We define the average number of local connections received by a neuron as K, 
and the external connections as -ft'ext- We assume that the connectivity matrices 
are random and we consider two scenarios (represented schematically in Figj2K,b): 

1) The network is fully connected {K = N, K^^t = Next) with random connection 
strengths (All-to-All, Fig|2^), characterized by a Gaussian distribution. 

2) The network is sparse, only a fraction of connections exists {k = K/N, 
kext = Kf,xt/Next)i selected at random but of constant strength (Sparse, Figj2)3), 
characterized by a Bernouillian distribution. 

The randomness of connections makes the network akin to a "disordered" sys- 
tem, which is characterized by a random but fixed substrate. We adopt a single 
notation for either case, All-to-All or Sparse network, by defining the mean and 
variance of matrix elements and their scaling with K, N: 

(G.,) = -kg/^fK (AGJ) = AViV (2) 

{Gf) = KxtQextl {K^xt (AGtf) = XlJN^xt. (3) 

(angular brackets denote average over the matrix distribution, A indicates vari- 
ation around the mean). In the All-to-All network, k = 1 and K = N, and the 
parameters of the matrix distribution are g and A. In the Sparse network, con- 
nection strengths have only two possible values, either zero or —g/y/K, the two 
parameters are g and k, but for convenience of notation we also use the parameter 
X"^ = g'^i^ ~ k). Similar equations and parameters describe the external connec- 
tions in either case. The mean connection is negative for G (inhibitory) and 
positive for Gext (excitatory), since g and gext are positive. Theoretical analysis 
also considers the case of local excitatory and inhibitory populations (in which 
(G) is a generic rank-1 matrix, see Appendix 2). 

We assume that the external activity :s.ext{t) is a stochastic process un corre- 
lated in both space and time, i.e. a white noise characterized by mean xf^*(t) = 
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Figure 2: Scheme of the network and distribution of activity and correlations. We 
study two network architectures: AU-to-All connectivity with random strengths 
(a), and Sparse random connections of fixed strength (b). Connection strength is 
illustrated by the thickness of edges. The distribution of activity across neurons 
(c,d) and the distribution of correlations across neuron pairs (e,f) are Gaussian 
for both types of networks ((c,e) for the All-to-AU network, (d,f) for the Sparse 
network). The parameters used are: g = gext = 1, Xext = Axg^.^ = 1, Kext = K. 
For the All-to-All network: K = 1000, Xext = 1, A = 1/^2. For the sparse 
network: K = 880 and k^xt = k = 1/2, which correspond to A = Xext = 1/v^- 



Text and covariance Axf^*(t)Axj^*(t') = /S.x1xi5ij5{t — t') (overline denotes the 
average over different realizations of the noise, and 5 denotes either the discrete 
Kronecker or continuous Dirac function). Therefore, Eq.(Il]) corresponds to a 
Ornstein-Uhlenbeck stochastic process [T6] . 

Results 

We study neural activity and correlations among neurons in a heterogeneous 
neural circuit model. Local recurrent connections are dominated by inhibition, 
while external feed-forward projections are excitatory. We consider two types 
of circuits, All-to-All connectivity with random strengths (Figj2^) and Sparse 
random connections of fixed strengths (Fig|2]3). Results are displayed in a single 
notation for either case (see Methods). Fig|2l:;,d shows the distribution of activity 
across neurons, and Figj2fe,f shows the distribution of correlations across neuron 
pairs. The purpose of this work is to describe how the mean and variance of those 
distributions depend on the parameters of the neural circuit. 

Since the model is linear, the neural activity x can assume negative values, 
while the activity of real neurons is positive (firing rate); Those values are inter- 
preted as deviations from a steady state, attained by the network in absence of 
the external input, around which the neural dynamics is approximately linear. If 
we denote the steady activity as /o, the firing rate is equal to / = /o -|- ex, where 
we assume that / is positive and c is a positive constant. As long as the linear 
approximation is valid, the correlations observed in the model are insensitive to 
the nature of the steady state (i.e. to the values of /o and c). 

Due to the linearity of the model, all quantities of interests can be simply 
calculated; The novel contribution of this work is averaging those quantities over 
the randomness of the connectivity matrix. Because connections are heteroge- 
neous, different neurons have a different activity, and we compute the sample 
mean across neurons in order to obtain the spatial average. If the number of 
neurons A^ is large, this is independent on the specific realization of the connec- 
tivity, therefore we perform its average over the distribution of connections, and 
we obtain (see Eq. flT3|) in Appendix 1; Angular brackets denote averaging over 
the random connectivity, overline denotes temporal average) 



The numerator of this expression is equal to the mean excitatory input received 
by a neuron, Qexty/Kext '^exti while the denominator expresses the recurrent inhibi- 
tion, whose total post-synaptic strength is gyK. Therefore, the strong recurrent 
inhibition counterbalances the large excitatory input and determines a relatively 
low activity, regardless of the network size. Note that the number of local and 
external connections, K and Kext-, are both large, but they tend to balance in 
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Figure 3: Mean (a) and Standard Deviation (b) of activity across neurons. 
The mean is plotted against the number of connections K, and the SD is plotted 
against the heterogeneity A. Analytical results (lines) are obtained from Eq.(j4]) in 
panel a and Eq. ([5]) in panel b. Simulation results are shown for the All-to-AU (blue 
circles) and Sparse (red triangles) network. Empty symbols show simulations of 
the neural dynamics, Eq.(IT]); Filled symbols show numerical evaluation of Eq. Olip . 
The parameters used are, in both panels, g = g^xt = 1, Xext = ^'^Ixt = 1- Iii 
panel b, K^xt = K = 500. For the All-to-AU network: Ag^j = 1, and A = l/-\/2 
in panel a. For the Sparse network: k^xt = k = 1/2 in panel a (which correspond 
to A = Xext = l/v2), while k^xt = k is varied in panel b. 



the expression above. Figj3^ shows an example of mean activity as a function of 
the number of connections; The mean activity is rather insensitive to the number 
of connections, which are taken equal to the external ones in each simulation 
[K = Kext)- The analytical result, Eq.(jl]), agrees with numerical simulations of 
the network dynamics, in both the All-to-All and the Sparse networks. 

Different neurons have different connections and therefore different activity, 
and the extent to which the activity varies from neuron to neuron is determined 
by the spatial variance. We calculate this quantity by taking the sample variance 
across neurons and averaging over the random connectivity, and we obtain (see 
Eq. (TT6l) in Appendix 1) 



Ax 



T^2 




■^^[{xf \^ + xl^^Xl^^ (5) 



The spatial variance of neural activity increases with the network heterogeneity, 
expressed by A and Xext for, respectively, the recurrent and external connections. 
Increasing the heterogeneity of connections increases the differences in the total 
input between neurons and therefore in their activities. The spatial variance is 
also proportional to the mean activity, local (x) and external Text- Furthermore, 
increasing the heterogeneity of recurrent connections leads to a divergence of 
the spatial variance, when A^ approaches one. In this case the state x = 0, 
which we assumed to correspond to a steady state of firing rates, destabilizes, 
and the linear approximation fails (see Appendix 1). FiglSb shows an example of 
the spatial variance as a function of the variability of the recurrent connections. 
The analytical result, Eq.(|5]), agrees with numerical simulations of the network 
dynamics, in both the All-to-All and the Sparse network. 

After looking at the mean and spatial variability of neural activity, we turn 
to the main theme of our work, the analysis of temporal variability and corre- 
lations. The activity of each neuron fluctuates in time, due to the fluctuating 
input, and those temporal fluctuations may be correlated since different neurons 
receive shared input. We study temporal variability and correlated fluctuations 
by calculating the covariance matrix, in particular the instantaneous covariance, 
at zero time lag. First, we look at the on-diagonal elements of this matrix, which 
are the temporal variances of different neurons. To determine the average tem- 
poral variance, we take the sample mean across neurons and we average over the 
random connectivity, obtaining (see Eq.( !T9|) in Appendix 1) 



(A^^ = (Ay-Axn = ^?2Ei 




'^extQext I Ag^j 



+ 



g 



Vk VT^T' 



(6) 



The temporal variance of neural activity is the sum of two pieces: the first term 
decreases with the number of connections as K'^^"^, while the second term remains 
finite. The first term indicates that recurrent inhibition (g) reduces temporal 
fluctuations. In fact, the inhibitory feedback not only reduces the mean activity 
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(Eq.(|4])), but also cuts down fluctuations by quickly counterbalancing the external 
excitatory input. This can be verified by calculating the instantaneous covariance 
between the external excitatory and the local inhibitory input, which is found 
large and negative, equal to —kextg1xt9^^- Furthermore, as in the case of spatial 
fiuctuations, temporal fiuctuations increase with the heterogeneity of connections, 
recurrent (A) and external {Xext}- Temporal fiuctuations also diverge when the 
network approaches the instability point, when the linear approximation fails 
(A^l). 

How much of the total variance, expressed in Eq.dH]), is independent rather 
than shared between neurons? To answer this question we calculate the average 
covariance, by looking at the off-diagonal elements of the covariance matrix, the 
pairwise covariances. We take the sample mean across neuron pairs and average 
over the random connectivity to obtain the average covariance (see Eq. (l20|) in 
Appendix 1) 







2 l + g^ 



Notably, this is equal to the first term in the total variance, Eq.([6]), implying 
that the two terms in the total variance express, respectively, the correlated and 
uncorrelated fiuctuations. Therefore, while the uncorrelated variance remains 
finite for large K, the correlated variance vanishes. The activities of neuron pairs 
tend to covary, due to their shared external input, but the recurrent inhibition 
makes the covariance small, of order K~^/'^ . In the Sparse network, the covariance 
vanishes if the probability of external connections is small {kext -^ 0), since the 
shared external input between neurons tends to zero in that case. In the All-to- 
AU network, the mean covariance vanishes if the mean input connection is zero 
{dext = 0); In that case, even if neurons receive a shared external input, neuron 
pairs may weight different inputs with the same or opposite signs, leading to 
respectively positive or negative covariance. Therefore, while the mean covariance 
across neuron pairs is zero, the covariance of single pairs may be positive or 
negative. 

The mean correlation is obtained by dividing the covariance, Eq.Q, by the 
variance, Eq.([6]), i.e. (we assume that variance and covariance are independent): 



Ax Ax ) 1 

(^) = ^7=^ = : I^^TTTF (8) 



Ax2) 1 



K-rt l+cl^/K 



ext 



This expression is positive and never exceeds one. It indicates that the mean 
correlation is small, of order K"^^"^, despite the strong and shared excitatory input 
between neurons. However, this result holds only in presence of the local recurrent 
inhibition [g > 0), and provided that external connections are heterogeneous 
(■^Li 7^ 0)- Heterogeneity of local connections (A) also contributes in decreasing 
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Figure 4: Mean correlation in the activity between neuron paris, Standand De- 
viation is shown in the insets. The correlation is plotted against the number of 
connections K (a) and the heterogeneity A (b) of the network. Analytical results 
(lines) are obtained from Eq.®. Simulation results are shown for the AU-to-All 
(blue circles) and Sparse (red triangles) network. Empty symbols show simula- 
tions of the neural dynamics, Eq.(IT]); Filled symbols show numerical evaluation of 
Eq.(TT7ll. The parameters used are, in both panels, g = g^xt = 1, Xext = ^xl^t = 1- 
In panel b, -R'ext = K = 500. For the All-to-AU network: Ae^i = 1, and A = l/-\/2 
in panel a. For the Sparse network: kext = k = 1/2 in panel a, which correspond 
to A = Xext = 1/ v2, while kext = k is varied in panel b. 
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the correlation. 

Therefore, the inhibitory feedback and the random connectivity are respon- 
sible for the small correlation. If the inhibitory feedback is removed, g = 0, the 
correlation becomes large. If the network heterogeneity is removed, A = Xext = 0, 
the correlation is equal to one, because the network is homogeneous and all neu- 
rons get the same input (but see Appendix 1 for a caveat on this case, Eq. fl22|) ). 
FigJH shows an example of the mean correlation as a function of the number of 
connections and the heterogeneity of the network. The analytical result, Eq.®, 
agrees with numerical simulations of the network dynamics, in both the AU-to-AU 
and the Sparse network. Insets in FigH] show the Standard Deviation of corre- 
lations, which appear to decrease with the number of connections as K~^^'^, and 
to increase with the heterogeneity of the network. 

For the Sparse network, it is interesting to note that all the above results 
depend on the number of neurons A^ only through the parameter A^ = (7^(1 — 
K/N). Therefore, while the above quantities may depend crucially the number 
of connections K, they depend only weakly on the number of neurons A^. In 
particular, the order of magnitude of correlations A'" 2 holds regardless of the 
number of neurons, which may be taken even infinite for any fixed value of K 
(provided that g < 1 and kext is non-zero). While in other studies the results 
are often described in terms on the number of neurons A^, here the determinant 
variable is the number of connections K. 

Discussion 

We found that inhibitory feedback and heterogeneous connections have important 
effects on the dynamics of the activity in a neural circuit. The strong excitatory 
input, shared between neurons, tends to drive the network to a highly active and 
correlated state. The inhibitory feedback is responsible for balancing the network 
activity, and also for reducing temporal fluctuations, in particular the correlated 
fluctuations across neurons. The heterogeneity of couphngs also plays a crucial 
role in reducing correlations, since homogeneous connections would determine 
homogeneous and therefore highly correlated activity. As a consequence, the 
observed mean correlation is positive and of small magnitude. The fact that mean 
correlation is positive is obvious, since neurons in a large population cannot be 
anti-correlated on averagqj. What is not obvious is that the mean correlation is 
of small magnitude. 

In presence of heterogeneous connections and inhibitory feedback, the mean 



^The mean correlation must be larger than —jfzi, therefore it must be positive in infinitely 
large populations. Proof: Any covariance matrix Q is positive definite, therefore h'^Qh > for 
any vector h. If we choose hi = 1/ ^/QTi and define the correlation matrix Rij ~ Qij / y^QuQjj 
we have that 
{R) > -jj^ 



we have that J2ij Rij = N + N{N — 1) (R) > 0. Therefore the mean correlation must be 
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correlation decreases with the number of connections as K~2. This order of mag- 
nitude agrees with experimental observations: a cortical neuron receives about 
K ~ 10^ connections [7] and the measured mean correlation is of order 10~^. 
Note that most experimental studies report a higher value of the mean corre- 
lation in the mammalian cortex, around 10~^ [H]. However, different studies 
measure correlations at different timescales, by integrating activity (or corre- 
lations) on different time windows. Larger temporal windows imply larger co- 
variance, because neurons covary on multiple timescales, and the increase in 
covariance is not accompanied by a comparable increase in the variance, due to 
after-hyperpolarization effects [1]. Therefore, a significant increase in the corre- 
lation, of about one order of magnitude, is observed when increasing the time 
window from about Ims to Is (see Fig{T]). Here we consider only instantaneous 
correlations: when measured at a short timescale (~ 1ms), the correlation is 
about 10~^. Besides the cerebral cortex, our model may be useful to investigate 
correlations in other brain areas, especially those dominated by inhibition, such 
as the striatum and globus pallidus. Interestingly, large correlations between pal- 
lidal neurons have been recently proposed as the origin of Parkinsonism [331 US]. 

Other modeling studies have addressed the issue of correlations in neural 
circuits. In the studies |35l |21] , the mean correlation was found to decrease with 
the number of neurons as A^~^ instead of the order K~2 that we found here. One 
difference between those models and ours is that we study a firing rate model 
instead of a spiking model, but that may not be a crucial difference. In those 
studies, the small correlation has been explained by an exact and instantaneous 
tracking of the excitatory input by recurrent inhibition. In particular, j35] points 
out that instantaneous tracking is the key for the iV~^ scaling of correlations, and 
if tracking is delayed then correlations would be of the order A^~2 instead. In our 
model, although recurrent inhibition is largely anti-correlated with the excitatory 
input, tracking is not instantaneous, since the excitatory input and the recurrent 
inhibition fluctuate on different timescales; Excitation is fast, while inhibition is 
slowed down by the recurrent dynamics. This difference may be responsible for 
the higher order of correlation observed in our study. 

While neural dynamics occurs on a variety of timescales in our modeo, as 
well as in real neurons |H], it may not be fast enough for instantaneous tracking. 
It would be instructive to slow down the fluctuations of the excitatory input in 
our model, for example by injecting colored instead of white noise in the net- 
work. In biological networks however, it is not known whether inhibition can 
track excitation exactly and instantaneously, as described in [351 121]. Strong 
anti- correlations between excitatory and inhibitory inputs have been observed 
experimentally [SHIEIEZ], but it remains unclear to what extent they may decor- 
relate neural activity. 



^The network is endowed with one very fast timescale, g ^K ^'^, while the remaining 
timescales are in between (1 + A)~^ and (1 — A)~^, in time units of r. 
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The difference in the order of magnitudes observed in our work compared to 
previous modehng studies may seem subtle, but this difference may be substantial 
from a functional standpoint. The magnitude of correlations found in our study 
is in between the "asynchronous" state [171 ES], in which the mean correlation 
is ~ A^~^, and the "synchronous" state, in which correlations are of order one. 
Those two models are opposite to one another, not only from a phenomenological 
point of view, but also in regard to the interpretation of the neural code and the 
hypotheses on the mechanisms underlying brain function. The former predicts 
that information is transmitted in the brain through firing rates, while the lat- 
ter suggests that is transmitted through synchronization of neural populations. 
There could be some truth in both views, but a satisfactory model bridging the 
two is not yet available. Further studies of correlations may provide the necessary 
understanding to fill the gap. 
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Appendix 1: Statistics of random networks 

In this section, we calculate the averages of neural activity and correlations with 
respect to both temporal fluctuations (noise) and the spatial variability of the 
connection strengths (disorder). Due to the linearity of the model, all quantities 
of interests can be simply calculated; The novel contribution of this work is 
averaging those quantities over the randomness of the connectivity matrix. The 
equation of dynamics ([1]) can be expressed in matrix form (the time constant of 
temporal evolution r is set to 1): 

^ = (G-/)x(t) + a..Xe,,(t) (9) 

where x is the vector of local neural activities, Xe^t is the vector of external 
neural activities, the matrices G (of size N x N) and Gext (size A^ x N^xt) express 
respectively the recurrent connections and the feed-forward projections, and / is 
the identity matrix. The equation of dynamics is linear and, given the interaction 
matrices G, Gext and the input signal Xe^t, the neural activity can be expressed 
as a sum over the external input weighted by an exponential temporal decay 



/t /'+00 

rft'e(^-^)(*-*')Ge.iXe.,(t') = / rft'e(^-^)*'a.txe.i(t - t') (10) 
-oo Jn 
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We assumed that initial conditions have decayed and that the inequahty A < 1 
holds, to prevent network activity from growing in time without boundqj. 

We start by calculating the mean neural activity. We perform the temporal 
average of the above expression, therefore we substitute the external activity 
Xext(^) with its average Text, and we perform the integral, obtaining (temporal 
average is denoted by overline) 

ll = Xext{I-GY^Gexa (11) 

where the vector 1 has all Next components equal to one. Because the matrices 
of connection strengths are heterogeneous, G and Gext-, different neurons have a 
different mean activity. In order to calculate the spatially averaged activity, we 
compute the sample mean across neurons. For large iV, this is independent on 
the specific realization of the spatial disorder, therefore we perform its average 
over the distribution of connectivity strengths, namely 

l?:) = (l-,y.^^ ^ (^'i-Hl ' G)-^G,^a) (12) 



The average (angular brackets) is across all possible realizations of the random 
matrices G and Gext- We denote by f the transpose operation. Note that, in 
the expression above, the row vector 1''' has A^ components, while the column 
vector 1 has Next- In the following we will use the same notation regardless 
of the dimension of 1, since that can be determined by the dimension of the 
multiplied matrix. Since G and Gext are independent, we can substitute Gext 
with its mean, {Gext) = Ss^i^^— s^nt. Furthermore, we show in Appendix 2, 

Eq.(|36]), that < !+(/ - G)-^l >= N{1 + gy/Ky\ Therefore, the mean activity 
is equal to 



; \ 9ext\ J^ext ( T o\ 

{x)-^^-^Xext (13) 

This expression is used in the main text (Eq.(j4])). 

Different neurons have different connections and therefore different activity, 
and the extent to which the activity varies from neuron to neuron is determined 
by the spatial variance. We calculate this quantity by taking the sample variance 
across neurons and averaging over the spatial disorder. We take the scalar product 
of Eq. Olip with itself, and we use again the fact that the sample mean does not 
depend on the spatial disorder for large A^, to obtain 



•^In the limit of large N, the real part of the eigenvalues of G is bounded by A. Therefore, 
if A > 1, some eigenvalues of (G — /) have non-negative real part, and the integral does not 
converge. This corresponds to an unstable fixed point at x = and network activity grows in 
time without bounds. 
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ltGL(/ - G^)-\I - G)-'Ge,tl ) - {xf (14) 



We also use the trace operator and its cyclic invariance, to rewrite this expression 
as 



Ax 



™2 




G't)-i(/-G)-^Ge..tll^GL 



{xf 



(15) 



Again, since G and Gext are independent, we can average separately the factors 
involving the two matrices. A simple calculation shows that (Gextll^Ge^.^ ) = 

(^g^ji^extll^ + ^Ixt^- Furthermore, we show in Appendix 2, Eqs.( l48|49|) that the 
following two equalities hold < Tr((/ - G)-\I - Gt)-i) >= n{1 - \^y\ and 
<Tr((/-G't)-i(/-G)-illt) >= N{l-\^y\l + gVK)~^. Using the expression 
of the mean activity, Eq. (ll3p . the spatial variance is equal to 



At 



I 



1-A2 



{xyy + xt,,x 



ext ^ ext 



(16) 



This expression is used in the main text (Eq.(l5])). 

After looking at the spatial variability, we study temporal variability and 
correlated fluctuations by calculating the covariance matrix. We take the scalar 
product of Eq.flTOl) with itself and we perform the temporal average, using the fact 
that the external stimulus is uncorrelated in space and time. This corresponds 
to the covariance matrix of a Ornstein-Uhlenbeck process |T6] , and is equal to 



Q = AxAxt 



Ax^ 



ext 



+ 00 



dt e (~j^ extXj' p-rt 6 



(17) 



Note that the covariance matrix satisfles the equation (G — I)Q + Q{G^ — /) = 
GextGexti but this cannot be used for averaging Q since G and Q are dependent 
and they do not commute. Also G and GextG\xt do not commute. The on- 
diagonal elements of the covariance matrix are the temporal variances of different 
neurons. To determine the average temporal variance, we take the sample mean 
across neurons and we average over the spatial disorder, obtaining 



N 



Ax' 



^EA-l 



Ax2 



ext 



i=l 



N 



[^ dt (Tr (e(«^"^)*e(^-^)*Ge.*GL) ) (18) 



where we applied the trace operator and used its cyclic invariance. A simple 
calculation gives (GextGext 
in Appendix 2, we obtain 



calculation gives (GextGext) = ^e^tS'Ltl-'-^ + '^Lt-^- Furthermore, using Eqs 05Of5ip 



Ax' 



Ax^ 



ext 



'^extQext 



2V1-A2[l + ^/^ 



ri + v/T^ 


A2(l + ^^ 


/K) 


y i + v^ 


-\^ + g^ 


fK 



+ ^ext 



(19) 
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In the main text, we substitute the term in square brackets with its hmit for large 
g^/K, which is equal to Vl — A^. The result is Eq.([6]). 

Next, we calculate the average covariance, by looking at the off-diagonal el- 
ements of the matrix in Eq.( lT71) . The off-diagonal elements are the pairwise 
covariances, and the sample mean across neuron pairs can be averaged over the 
spatial disorder to obtain the average covariance 



, / 1 1'^ \ 

AxAx) = ( — -VAxjAa;,- 






Using Eq.f l52|) in Appendix 2, we obtain 



^"^ / ^ dt (Tr (e(«^-^)*(llt - J)e(«-^)*a..GL 



A^A^\ = ^^ ^--tgl^t (20) 

This expression is used in the main text (Eq.(j7])). 

The mean correlation is obtained by dividing the covariance, Eq. fl20|) . by the 
variance, Eq.f llQp . i.e. (we assume that variance and covariance are independent): 



Ax Ax) J\ - \2 

/ Ax2\ 1+-/T^^{1+9VK) , Ag^,(l+gv^) '^ ^ 

This expression is positive and never exceeds one. In the main text, we substitute 
the first term in the denominator with its limit for large gvK, which is equal to 
Vl — A^. The result is Eq.(IH]). Note that, due to this simplification, the equation 
in the main text may be inaccurate in some cases. For example, if the external 
heterogeneity is removed, Xext = 0, then Eq.® gives {R) = 1, while Eq. (l2T]) gives 

{R) = l , ^ (22) 

This result is more accurate since, if the local heterogeneity is not removed, 
A 7^ 0, then the network is not precisely homogeneous and the correlation is not 
expected to be exactly equal to one (however, correlation is one for large K even 
in Eq. (l22p ). The two expressions also differ for g = 0. 

Appendix 2: Traces of random matrix products 

In this section we introduce the diagrammatic notation to calculate the quenched 
averages of random matrix products (see e.g. [2U]). Theoretical results are ob- 
tained for the Gaussian distribution, although numerical simulations suggest that 

16 



they generalize to other distributions with the same mean and variance (e.g. 
BernouiUi). We will consider the case in which the mean of the matrix element is 
~ g/N and then recover the scaling studied in the main text by analytical con- 
tinuation and the substitution g — )■ gyK. We conclude the section by studying 
the case of non-homogeneous mean (e.g. interconnected excitatory and inhibitory 
neurons) . 

We start with the problem of calculating the quenched average of the trace 
of a power of the random matrix R in the limit of large A^ (where the size of the 
matrix is iV x N). The matrix R is characterized by independent and normally 
distributed elements, each element having zero mean and variance A^~^, namely 

,2 V 1 



(i?,,) = (4) = - (23) 

We start by calculating the second order, i.e. the average trace of the square of 
R. For convenience of notation, we omit the sum over the indices (in this case 
the sum over the indices a, 6, c, d) 

Tr (i?2) \ = dad^bc (RabRcd) = N-' SadSbJacSM (24) 



The diagram corresponding to this expression is shown in FigJS^. The diagram 
is obtained by drawing one node for each one of the four indices a, b, c, d, and 
by drawing an edge for each delta function in the expression, where the two 
nodes connected by the edge correspond to the two indices of the delta func- 
tion. Horizontal edges are due to the operations of trace (base edge) and matrix 
multiplication (middle edge), while arc-shaped edges are due to averaging. The 
multiple edges determine different paths, and each pair of nodes connected by a 
path (even if not linked by an edge), corresponds to a pair of indices that must 
be equal, since they are connected by a sequence of delta functions. Therefore, 
for each closed loop in the diagram there is one redundant delta function, which 
can be eliminated without performing the sum over the corresponding indices. 
This implies that each closed loop contributes with a factor N, due to a free sum 
over A^ elements. Since the diagram for the second order has one loop, we have 
(Tr(i?2)) = N-'^N = 1. 

Note that all terms of odd order are zero, because (Rij) = for odd k. The 
next order is therefore the fourth order, which is equal to (again we omit the sum 
over all indices) 

Tr (i? ) ) = SahSbc^de^fg (RabRcdRefRgh) = (25) 



= N 6ahSbcSdeSfg [SacSbdSegSfh + ^aehf^cg^dh + ^aghh^ce^df] (26) 

The fourth order has three diagrams, one for each term in the sum, shown in 
FigJSb- The middle diagram has two closed loops, while the other two have 
only one loop. Therefore the other two terms can be neglected, the middle term 
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Figure 5: Diagrams of the traces of random matrix powers, described by Eq. fl27|) . 
Below each diagram the number of its closed loops is indicated, a) Second order, 
b) Fourth order, c) Sixth order. 

contributes with a factor A^^ and the fourth order gives (Tr(i?^)) = 1. The 
contribution of the fourth order moment (( i?f- )) can be neglected. 

The fifteen diagrams of the sixth order are shown in FigJSt. Again, we neglect 
moments higher than the second, and we note that only one diagram contributes 
with three loops, therefore (Tr(i?^)) = 1. By iterating this procedure, we find 
that order 2k has {2k — 1)!! diagrams of which only one has k loops, therefore 

)2fc 



Tr(i?^'^j) = l (27) 

for all values of k. 

Note that the elements of the matrix R have zero mean, while the matrices 
considered in the main text {G and Gext) have non-zero mean. As explained 
below, in order to calculate the average trace of matrix powers with non-zero 
mean, we need to compute averages where R is interleaved by the matrix of ones. 
We denote by 1 the column vector of A^ components all equal to one, by 1^ the 
row vector, and by 11^ the N x N matrix with all elements equal to one (we 
denote by f the transpose operation). We consider the two second order terms 



Tr [Rll^R) ) = 6ad (RabRcd) = N-^5ad5aAd = 1 (28) 

Tr (i?2llt) \ = 3,, (RabRcd) = N-'6f,cSaJM = 1 (29) 



It is not surprising that these two expressions are equal, since the trace is cyclic 
invariant. The only difference of these expressions with Eq. fl2^ is the absence of 
a factor 6bc in the former expression, and 6ad in the latter. This corresponds to 
cutting, respectively, the middle and the base horizontal edges in the diagram of 
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FiglS^. In general, inserting a matrix of ones at a given point of the sequence 
of R products is equivalent to cutting the horizontal edge at that point in the 
corresponding diagram. If the edge belongs to a closed loop, the cut has the 
only effect of removing a redundant delta function, and there is no change in the 
contribution of that diagram to the sum; Conversely, if the edge belongs to an 
open path, the cut determines an additional A^ factor, because the delta function 
removed was not redundant. Since all diagrams have at least one closed loop, 
inserting a single matrix of ones has no effect at all orders. Therefore, 



Tr(i?2 



ll^R''')) = l (30) 



for all k' = 0, . . . , 2k. Unless more loops are available to cut, inserting more 
matrices of ones may cut open paths, therefore the trace may be multiplied by 
A^. An additional A^ factor is obtained also by multiplying the matrix of ones 
with itself, which occurs whenever additional matrices are inserted at same point 
in the sequence (we have that l'''! = A^ and (11''')^ = A^'^"^!!''' if A; > 0). 

Using the above results, we can calculate the average trace of random matrix 
powers with non-zero mean and arbitrary variance (provided that the variance is 
of order N~^). We consider the matrix G equal to 

G = j-ll^ + XR (31) 

Note that the mean of this matrix has a different scaling with respect to that 
considered in the main text, but we will recover the latter by the substitution 
g — >■ —y/Kg. A power of G is calculated by multiplying G to itself, and this 
determines an ordered product of powers of the matrices R and 11^. Note that 
these two matrices do not commute, therefore the binomial theorem cannot be 
applied. We consider the average trace 

(lY {g') ) = Y: N-^'g^'X'^'' E (Tr(. . .)) (32) 

fc'=0 

where the trace in the right hand side is applied to an ordered product of k' 
matrices 11^ and k — k' matrices i?, and the sum runs over all the L ,) ordered 
products for a given k and k' . Using the above results, we find that the contri- 
bution of any of those traces is zero for k — k' odd, is equal to one for k' = 
(provided that k is even), is equal to A^'^ for k' = k and is at most of order A^^ ~^ 
for k' = 1, . . . ,k — 1. Therefore, the leading order terms are k' = k (for any 
value of k) and k' = (for k even), and all other terms can be neglected. If the 
matrix G'^ is further multiplied by a matrix of ones, the term k' = can also be 
neglected, and we find that 

Tr (g'^II^) ) = Ng^ = Tr ({Gf 11^) (33) 
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Figure 6: Diagrams of the traces of random matrix powers multiplied by its 
transpose, described by Eq. (l40|) . Below each diagram the number of its closed 
loops is indicated, a) Second order, b) Fourth order, c) Sixth order. 

for all values of k. Note that if the mean of G has a higher order in N, the 
result still holds. This expression is particularly useful to compute the average 
of bracket expressions. Because Tr (Axy'^]=y'^Ax for any matrix A and vectors 
X, y, the expression can be rewritten as 

(i^GH) = 1^G)^1 (34) 

for all values of k. Since any infinitely different iable function / can be expanded 
in Taylor series, the above result implies that 

lt/(G)l) = lt/((G))l (35) 



Therefore, the following expression can be calculated and used to compute the 
mean activity in the main text. 

It (/ - Gy' l) = T^ (36) 

/ 1- g 



Note that the substitution g -^ —gvK must be applied to recover the scaling 
studied in the main text. 

Next, we calculate the diagrammatic expansion for products of a random 
matrix with its transpose. Again, all odd orders vanish and we neglect moments 
higher than the second at all orders. The second order term is 

Tr (RR^) ) = 6adSbc (RabRdc) = N'' SadSbJadSbc (37) 
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The corresponding diagram has two loops and is shown in FigjG^, therefore the 
loops contribute with a factor A^^ and the second order is ( Tr (RR^j ) = N. The 
fourth order is equal to 



Sah^bc^deSfg {RabRcdRfeRhg) — (3^ 



Tr (r^R 

= N^"^ SahSbcSdeSfg [SacSbdSfhSeg + ^afhe^ch^dg + ^ah^bg^cf^de] (39) 

The three diagrams are shown in Figl6]D. The first two diagrams have one loop, 
while the third has three. Therefore, that diagram contributes with a factor N^ 
and the fourth order is equal to (Tr (R^Rj ) = N. The diagrams for the sixth 
order are shown in FiglHb: only one diagram has four loops, and no diagram has 
three, therefore (^Tr (R^R'^ j) = N. Iterating the procedure, we find that order 
2k has {2k — 1)!! diagrams of which only one has k + 1 loops, therefore 

Tr (r'^R''^) ) = N (40) 



for all values of k. Other combinations of powers of R and its transpose give a 
smaller contribution, i.e. ( Tr f i?^'^"'^ R'' )) ~ '^(^) ^°^ ^' t^ ^• 

Inserting matrices of ones in this case has a similar effect as in the case above, 
Eq. fl30l) . each matrix cuts the horizontal edge corresponding to where the matrix 
is placed. Again, since each diagram has at least one loop, the insertion of a single 
matrix of ones (and the consequent edge removal) has no effect on the trace at 
all orders. Therefore, 

Tr (r''-''' 11^ R''' R''^)) = (ti (R'^R'^'hl^R''-'''^)) = N (41) 



An insertion in a term with unequal powers of R and R^ remains of order one. 
Adding more matrices increases the trace by an order A^ for each matrix, provided 
that no further loops are cutted. 

Using the above expressions, we can compute the average of products of pow- 
ers of the matrix G and its transpose, namely 

k'=Ol'=0 

where the trace in the right hand side is applied to an ordered product of k' + /' 
matrices 11^, k — k' matrices R and / — /' matrices R^. li k' = k and /' = /, the 
trace is equal to A^'^'"'"', and the term is of order one. li k' = and /' = 0, the trace 
contributes with an order A^, provided that k = I. li k — k' = I — I', the traces 
contribute at most with an order A^'^ "*"' , and the term is of order one, while if 
k — k'^l — l' the term is of smaller order. Therefore the leading order is A^, and 
we have 
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[Tt (^G'^G^^)) = NdkiX''^^ (43) 

In the case in which matrices of ones are inserted, the term k' = 0, I' = is 
no longer leading, and many other terms have to be considered. Those are the 
terms for k — k' = I — I', and for which additional inserted matrices cuts the same 
loop. Since the leading diagrams at all orders have one two-nodes loop in the 
middle and one at the boundaries, if a matrix of ones is inserted in the middle or 
at the boundaries, additional matrices must continue to be inserted at the same 
place in order to cut the same loop. We eliminate one sum and we use the index 
m = k — k' = I — I' in place of k' and /'; We obtain 

min(fc,Z) 
m=0 

(44) 

mill (A:,?) 
(Tr (G^II^G^^) ) = J2 iV2m-fc-i^i+fc-2m,^2m ^r^^ ^-^n,^^^f^k+l-2m+l j^m.^^^ ^ 

m=0 

(45) 
Both expressions are equal to 

inin {k,l) 

(Tr (G^G'^II^)) = (Tr (G^II^G^^)) = N J2 d^^-'^'^X'-'^ (46) 

m=0 

Furthermore, we calculate the average trace with two inserted matrices. In that 
case, the leading term is for k = k' and / = /' (or m = 0), and we obtain 

Tr (GHl^G^hl^)) = N"" /+' = Tr ({G)'' llt(G)'^ll^') (47) 

Using the expressions above and the Taylor series expansion of infinitely dif- 
ferentiable functions, we calculate the following traces that are used in the main 
text to compute the variance and covariance of the activity 

Tr((/-G)-i(/-Gt)-i)) = ^ (48) 

Tr ((/ - G)-llt(/ - Gt)-) ) = (^_,,^^_^^, (49) 

'rfte--(Tr(e^*e^^*)) = ^-^ (50) 



dt e"2* (Tr (e^^llte^^*)) = , ^, 

\ ^ J/ 2VT^X^{l-g) 
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i + v^r^A^- 



(51) 



/;*e-(Tr(e«llte="llt)) = 5J^ (52) 

Note that the substitution g — )■ —g\fK must be apphed to recover the scahng 
studied in the main text. If K is proportional to iV, this substitution may change 
the order of magnitude of various terms in the summation considered above, 
possibly modifying the leading terms in each sum. Note that all series converge 
only for \g\ < 1, but their sum can be evaluated at (7 — ?■ —g\fK by analytical 
continuation. Then, approximating the sums by the leading terms described 
above is accurate under the assumption that all series involving lower order terms 
converge to bounded functions of g. 

We conclude this section by noting that the above calculations can be per- 
formed also in the case of a non- homogeneous mean. We have assumed that the 
mean matrix is homogeneous, i.e. is proportional to the matrix of ones (G) oc 11^. 
However, the same approach could be taken to analyze the more general case in 
which the mean is a rank-1 matrix, i.e. (G) = ab"!", where a and b are two arbi- 
trary vectors. For example, if the mean connection strength only depends on the 
pre-synaptic neuron, (G) = lb"'', then interconnected excitatoty and inhibitory 
neurons can be considered, where the vector b is positive for excitatory neurons 
and negative for the inhibitory ones. In that case, all results above still hold, 
with the simple substitution g — )■ A^~^b"''l. 
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